(Notebook version: 1.0 Initial - 2021 05 08)
The main goal of this notebook is to implement a simple, easy to use simulation tool to simulate FM signals for passive radar purposes. In the current state of the notebook and the provided API the generated FM signal is not compatible with the "ITU BS.450: Transmission standards for FM sound broadcasting at VHF" [1] recommendation. The implementation of this standard is a later task of this work.
In the most simple case an FM modulated broadcast signal can be generated as follows:
$$s(t) = sin \left(2 \pi f_c t - 2 \pi k_{fm} \int_0^t{s_{mod}(\tau)} d\tau \right),$$where $f_c$ is the carrier frequency and $s_{mod}(t)$ is modulating signal. The frequency deviation is defined as $f_d = k_{fm} max\{s_m(t) \}$, which is $75 kHz$ in the standardized FM broadcast stations.
To prepare the simulated data digitally we also need the discrete time version of the above defined signal. $$s[n] = sin \left(2 \pi \frac{f_c}{f_s} n - 2 \pi \frac{k_{fm}}{f_s} \sum_{m=0}^n{s_{mod}[m]} \right),$$
In the few next code section let us generates an FM signal with a symetrical square wave modulating signal.
# Python package imports
# For math calculations
import math
import numpy as np
# For plotting
import plotly.graph_objects as go
from plotly.offline import init_notebook_mode
init_notebook_mode()
#Default figure sizes
fig_w = 900
fig_h = 750
# Generate the modulation signal
# Simulation parameters
T_sim = 0.1 # simulation time [s]
fs = 100*10**3 # sampling frequency of the simulation [Hz]
Ts = 1/(5*10**3) # Square wave period
# Calculation
N = int(T_sim/(1/fs)) # Number of samples to simulate
s_mod = []
pol = -1
for i in range(math.ceil(N/(Ts*fs))):
s_mod.append(pol)
pol *= -1
s_mod = np.kron(np.array(s_mod), np.ones(int(round(Ts*fs)))) # Expand the signal to fit the timescale
s_mod = s_mod[0:N] # Crop extra signal samples
# Plot the modulating signal
print("Total number of simulated samples : {:d}".format(N))
print("Symbol period : {:.2f} ms".format(Ts*1000))
n = 100# number of samples to plot
fig=go.Figure()
fig.add_trace(go.Scatter(x=np.arange(n)*(1/fs)*1000,
y=s_mod[0:n],
mode='lines',
name=("Modulating signal")))
fig.update_layout(
autosize=False,
width=fig_w,
height=fig_h,
yaxis=dict(title_text="Amplitude "),
xaxis=dict(title_text="Time [ms]"),
title={'text': 'Symmetrical square wave modulating signal',
'y':0.9,
'x':0.5,
'xanchor': 'center',
'yanchor': 'top'}
)
Total number of simulated samples : 10000 Symbol period : 0.20 ms
# Generate the FM signal
fc = 10 # carrier frequency [kHz]
fd = 5 # frequency deviation [kHz]
k_fm = fd*10**3/np.max(s_mod)
t = np.arange(np.size(s_mod)) # sample indices
s = np.sin(2*np.pi*fc*10**3/fs*t - 2*np.pi*k_fm/fs*np.cumsum(s_mod))
# Plot the time domain modulated signal
n = 100# number of samples to plot
fig=go.Figure()
fig.add_trace(go.Scatter(x=np.arange(n)*(1/fs)*1000,
y=s[0:n],
mode='lines',
name=("Modulated signal")))
fig.add_trace(go.Scatter(x=np.arange(n)*(1/fs)*1000,
y=s_mod[0:n],
mode='lines',
name=("Modulating signal")))
fig.update_layout(
autosize=False,
width=fig_w,
height=fig_h,
yaxis=dict(title_text="Amplitude "),
xaxis=dict(title_text="Time [ms]"),
title={'text': 'FM signal, symmetrical square wave',
'y':0.9,
'x':0.5,
'xanchor': 'center',
'yanchor': 'top'}
)
Let us inspect the spectrum of the generated FM signal.
sw = np.fft.fft(s)
sw_abs = np.abs(sw)
sw_log = 20*np.log10(sw_abs)
sw_norm = sw_log - np.max(sw_log)
sw_align = np.fft.fftshift(sw_norm)
freqs = np.fft.fftshift(np.fft.fftfreq(N, 1/fs)) /10**3
fig=go.Figure()
fig.add_trace(go.Scatter(x=freqs,
y=sw_align,
mode='lines',
name=("Modulated signal")))
fig.update_layout(
autosize=False,
width=fig_w,
height=fig_h,
yaxis=dict(title_text="Amplitude [dB]"),
yaxis_range=[-50,0],
xaxis=dict(title_text="Frequency [kHz]"),
title={'text': 'FM signal, symmetrical square wave',
'y':0.9,
'x':0.5,
'xanchor': 'center',
'yanchor': 'top'}
)
In the following code sections let us generate a real world FM signal. For this we are going to use an imported music file whhich will be the modulating signal. For the algorithm development the first sound file will be a simple 1kHz sinusoid.
# Importing the sound file samples
from pydub import AudioSegment
from pydub.utils import mediainfo
fname = "music_samples/2k_test.mp3"
sound = AudioSegment.from_mp3(fname)
sound_samples = np.array(sound.get_array_of_samples())
info = mediainfo(fname)
fs_sound = int(info['sample_rate'])
print("{:d} samples have benn imported from the sound file".format(np.size(sound_samples)))
print("Sampling rate: {:d} Hz".format(fs_sound))
2654208 samples have benn imported from the sound file Sampling rate: 44100 Hz
Let us check that the importing is working properly.
# Plot the time domain signal
n = 5000 # Number of samples to plot
fig=go.Figure()
fig.add_trace(go.Scatter(x=np.arange(n)*(1/fs_sound)*1000,
y=sound_samples[0:n],
mode='lines'))
fig.update_layout(
autosize=False,
width=fig_w,
height=fig_h,
yaxis=dict(title_text="Amplitude "),
xaxis=dict(title_text="Time [ms]"),
title={'text': 'Imported 2KHz, single tone signal',
'y':0.9,
'x':0.5,
'xanchor': 'center',
'yanchor': 'top'}
)
To create the simulation we are going to cut out a signal portion from the full imported signal array. This cropped signal will has the same length as the simulation has.
At this point we have to take into consideration the required bandwidth of the simulated signal and the sampling frequencies as well. The FM broadcast signal has a frequecy deviation of $f_d = 75 kHZ$. According to Carson's formula: $$ B = 2\left(f_D+max\{f_m\}\right),$$ where $B$ is the FM signal bandwidth and $f_m$ is the modulating frequency. Considering that the maximum useful frequency of a sound wave is less than $20 kHz$, $f_m \leq 20 kHz$. Thus, $$B \leq 200 kHz$$
According to the Nyquist sampling theory the sampling frequency of the modulated signal should be at least $400 kHz$. $f_{s,FM} \geq 2B$. In order to work with a common sample rate when simulating the FM signal, we have to upsample the imported signal. $$ P f_{s, imp} = f_{s,FM},$$ where $K$ is the upsampling ratio, $f_{s, imp}$ is the sampling rate of the imported modulating signal.
import resample
# Parameters
P = 5 # upsample ratio
T_sim = 0.1 # simulation time [s]
offset = 10000 # sample offset
N_raw = math.ceil(T_sim/(1/fs_sound)) # Number of samples needed from the modulating signal
# Upsample the imported signal
sound_samples = sound_samples[offset:offset+N_raw] # Cut out useful portion
s_mod = resample.resample_n_filt(P,1,10,sound_samples)
N = np.size(s_mod)
fs_sound = P * fs_sound
s_mod = s_mod/np.max(np.abs(s_mod)) # Normalize
print("Number of samples to simulate {:d}".format(N))
Number of samples to simulate 22050
Let us generate now the baseband modulated signal
fs = fs_sound
fc = 0 # carrier frequency [kHz] -> Baseband
fd = 75 # frequency deviation [kHz] -> Same as in FM broadcast
k_fm = fd*10**3/np.max(s_mod)
t = np.arange(np.size(s_mod)) # sample indices
s = np.sin(2*np.pi*fc*10**3/fs*t - 2*np.pi*k_fm/fs*np.cumsum(s_mod))
# Plot the time domain modulated signal
n = 1000 # number of samples to plot
fig=go.Figure()
fig.add_trace(go.Scatter(x=np.arange(n)*(1/fs)*1000,
y=s[0:n].real,
mode='lines',
name=("Modulated signal")))
fig.add_trace(go.Scatter(x=np.arange(n)*(1/fs)*1000,
y=s_mod[0:n].real,
mode='lines',
name=("Modulating signal")))
fig.update_layout(
autosize=False,
width=fig_w,
height=fig_h,
yaxis=dict(title_text="Amplitude "),
xaxis=dict(title_text="Time [ms]"),
title={'text': 'FM signal, imported sine wave',
'y':0.9,
'x':0.5,
'xanchor': 'center',
'yanchor': 'top'}
)
sw = np.fft.fft(s)
sw_abs = np.abs(sw)
sw_log = 20*np.log10(sw_abs)
sw_norm = sw_log - np.max(sw_log)
sw_align = np.fft.fftshift(sw_norm)
freqs = np.fft.fftshift(np.fft.fftfreq(N, 1/fs)) /10**3
fig=go.Figure()
fig.add_trace(go.Scatter(x=freqs,
y=sw_align,
mode='lines',
name=("Modulated signal")))
fig.update_layout(
autosize=False,
width=fig_w,
height=fig_h,
yaxis=dict(title_text="Amplitude [dB]"),
yaxis_range=[-50,0],
xaxis=dict(title_text="Frequency [kHz]"),
title={'text': 'FM signal spectrum, sine wave',
'y':0.9,
'x':0.5,
'xanchor': 'center',
'yanchor': 'top'}
)
Let us now inspect the monostatic ambiguity function of the generated FM signal. For this, we are going to use the pyapril library.
from pyapril import detector
from pyapril import RDTools
max_speed = 1000 # maximum monostatic target speed [km/h]
max_range = 100 # maximum bistatic target range [km]
frf = 100 * 10**6 # radio frequency of the carrier [Hz]
print("Sample, number",N)
wavelength = 3*10**8/frf
max_Doppler = int(2* max_speed /3.6 / wavelength)
max_range_cell = int((max_range*1000)/(3*10**8)*fs)
chi = detector.cc_detector_fd(s, s, fs, max_Doppler, max_range_cell, verbose=0, Qt_obj=None)
# Plot the ambiguity function
fig = RDTools.plot_rd_matrix(chi, max_Doppler=max_Doppler, fs=fs)
fig.show()
Sample, number 22050
The following cell wraps the above developed scripts into a single function, so we can generate and analyze FM signals in a more convinient way.
def generate_fm_signal(s_mod, fs):
"""
Generates FM signal using the give modulation signal
Parameters:
-----------
:param: s_mod: Modulating signal
:param: fs: Sampling frequency of the modulating signal [Hz]
:type: s_mod: One dimensional numpy array
:type: fs: float
Retrun values:
--------------
:return: s : Generated FM signal
:return: fs: Sampling frequency of the modulated signal
:rtype: s: One dimensional numpy array
:rtype: fs: float
"""
# Internal processing parameters
fd = 75 # frequency deviation [kHz] -> Same as in FM broadcast
# Generate FM signal
s_mod = s_mod/np.max(np.abs(s_mod)) # Normalize
k_fm = fd*10**3/np.max(s_mod)
s = np.sin(2*np.pi*k_fm/fs*np.cumsum(s_mod)) # Modulate
return s, fs
def generate_fm_signal_from_sound(sound_fname, T_sim, offset):
"""
Generated FM signal from sound file
The offset parameter is a float number between 0 and 1. It
specifyies the start position of the time window(T_sim)
that will be used from the sound file for the FM signal
preparation.
Parameters:
----------
:param: sound_fname: Name of sound sound file to be imported
:param: T_sim: Duration of the reqested simulated signal [s]
:param: offset: Offset position of the processed window
:type: sound_fname: string
:type: T_sim : float
:type: offset: float (0..1)
Retrun values:
--------------
:return: s : Generated FM signal
:return: fs: Sampling frequency of the modulated signal
:rtype: s: One dimensional numpy array
:rtype: fs: float
"""
# Internal processing parameters
resample_ratio = 5
fir_tap_size = 10 # Resample filter tap size
fd = 75 # frequency deviation [kHz] -> Same as in FM broadcast
# Import sound file
sound = AudioSegment.from_mp3(sound_fname)
sound_samples = np.array(sound.get_array_of_samples())
info = mediainfo(sound_fname)
fs = int(info['sample_rate'])
# Resample
offset = int(offset*np.size(sound_samples))
N_raw = math.ceil(T_sim/(1/fs)) # Number of samples needed from the modulating signal
sound_samples = sound_samples[offset:offset+N_raw] # Cut out useful portion
s_mod = resample.resample_n_filt(resample_ratio,1,fir_tap_size,sound_samples)
fs = resample_ratio * fs
# Generate FM signal
s_mod = s_mod/np.max(np.abs(s_mod)) # Normalize
k_fm = fd*10**3/np.max(s_mod)
s = np.sin(2*np.pi*k_fm/fs*np.cumsum(s_mod)) # Modulate
return s, fs
def plot_fm_ambiguity(s, fs, frf, max_speed = 1000, max_range=100, **kwargs):
"""
:param: s: Illuminator signal
:param: fs: Sampling frequncy of the input signal
:param: frf: Carrier frequency of the FM signal [MHz]
:param: max_speed: [km/h]
:param: max_range: [km/h]
"""
dyn_range = kwargs.get('dyn_range')
wavelength = 3*10**8/(frf*10**6)
max_Doppler = int(2* max_speed /3.6 / wavelength)
max_range_cell = int((max_range*1000)/(3*10**8)*fs)
chi = detector.cc_detector_fd(s, s, fs, max_Doppler, max_range_cell, verbose=0, Qt_obj=None)
# Plot the ambiguity function
fig = RDTools.plot_rd_matrix(chi, max_Doppler=max_Doppler, fs=fs, dyn_range=dyn_range)
return fig
fm_sig, fs = generate_fm_signal_from_sound('music_samples/2k_test.mp3', T_sim=0.1, offset=0.01)
fig = plot_fm_ambiguity(fm_sig, fs=fs, frf=100, max_speed=1000, max_range=100)
fig.update_layout(
autosize=False,
width=fig_w,
height=fig_h,
title={'text': 'FM signal Ambiguity Function, 1 kHz sine wave',
'y':0.9,
'x':0.5,
'xanchor': 'center',
'yanchor': 'top'}
)
In this section we are going to calculate and display the ambiguity function of some typical FM signal, such as:
The total simulation time for all the diferent modulating singal types will be $T_{sim} = 1s$. Accordingly the maximum achvieable coherent processing gain is $G_{proc} \leq BT \simeq 200 kHz ~ 1s \rightarrow 53 dB$. Based on this rough estimation the dynamic range of the ambiguity functions plots are constrainted to $55 dB$ to highlight the relevant features of the signals.
T_sim = 0.5# [s]
plot_dyn_range = 55 # [dB]
# Parameters
fs = 200 *10**3
f0 = 10000
# Sine wave
s_mod = np.sin(2*np.pi*f0/fs*np.arange(int(T_sim*fs)))
# FM signal generation
fm_sig, fs = generate_fm_signal(s_mod, fs)
# Ambigutity function plot
fig = plot_fm_ambiguity(fm_sig, fs=fs, frf=100, max_speed=1000, max_range=100, dyn_range=plot_dyn_range)
fig.update_layout(
autosize=False,
width=fig_w,
height=fig_h,
title={'text': 'FM signal Ambiguity Function, {:.0f} kHz simulated sine wave'.format(f0/1000),
'y':0.9,
'x':0.5,
'xanchor': 'center',
'yanchor': 'top'}
)
# Parameters
fs = 200 *10**3
Ts = 1/(5*10**3) # Symbol period
# Square wave
N = int(T_sim/(1/fs))
s_mod = []
pol = -1
for i in range(math.ceil(N/(Ts*fs))):
s_mod.append(pol)
pol *= -1
s_mod = np.kron(np.array(s_mod), np.ones(int(round(Ts*fs))))
s_mod = s_mod[0:N]
# FM signal generation
fm_sig, fs = generate_fm_signal(s_mod, fs)
# Ambigutity function plot
fig = plot_fm_ambiguity(fm_sig, fs=fs, frf=100, max_speed=1000, max_range=100, dyn_range=plot_dyn_range)
fig.update_layout(
autosize=False,
width=fig_w,
height=fig_h,
title={'text': 'FM signal Ambiguity Function, {:.0f} kHz simulated symmetrical square wave'.format(5),
'y':0.9,
'x':0.5,
'xanchor': 'center',
'yanchor': 'top'}
)
# Parameters
fs = 200 *10**3
Ts = 1/(5*10**3) # Symbol period
# FSK data
s_mod=np.round(np.random.uniform(0,1,math.ceil(N/(Ts*fs))))*2-1
print("FSK data (sample): {0}".format(s_mod[0:10]))
s_mod = np.kron(np.array(s_mod), np.ones(int(round(Ts*fs))))
s_mod = s_mod[0:N]
# FM signal generation
fm_sig, fs = generate_fm_signal(s_mod, fs)
# Ambigutity function plot
fig = plot_fm_ambiguity(fm_sig, fs=fs, frf=100, max_speed=1000, max_range=100, dyn_range=plot_dyn_range)
fig.update_layout(
autosize=False,
width=fig_w,
height=fig_h,
title={'text': 'FM signal Ambiguity Function, {:.0f} kHz simulated FSK signal'.format(5),
'y':0.9,
'x':0.5,
'xanchor': 'center',
'yanchor': 'top'}
)
FSK data (sample): [ 1. 1. -1. 1. -1. -1. 1. -1. -1. 1.]
# Parameters
fs = 200 *10**3
# Gaussian noise
N = int(T_sim/(1/fs))
s_mod = np.random.normal(0,1,N)
# FM signal generation
fm_sig, fs = generate_fm_signal(s_mod, fs)
# Ambigutity function plot
fig = plot_fm_ambiguity(fm_sig, fs=fs, frf=100, max_speed=1000, max_range=100, dyn_range=plot_dyn_range)
fig.update_layout(
autosize=False,
width=fig_w,
height=fig_h,
title={'text': 'FM signal Ambiguity Function, Gaussian noise',
'y':0.9,
'x':0.5,
'xanchor': 'center',
'yanchor': 'top'}
)
# Generate FM signal
fm_sig, fs = generate_fm_signal_from_sound('music_samples/Für Elise (Piano version).mp3', T_sim=T_sim, offset=0.3)
# Ambigutity function plot
fig = plot_fm_ambiguity(fm_sig, fs=fs, frf=100, max_speed=1000, max_range=100, dyn_range=plot_dyn_range)
fig.update_layout(
autosize=False,
width=fig_w,
height=fig_h,
title={'text': 'FM signal Ambiguity Function, Für Elise (Piano version)',
'y':0.9,
'x':0.5,
'xanchor': 'center',
'yanchor': 'top'}
)
# Generate FM signal
fm_sig, fs = generate_fm_signal_from_sound('music_samples/ACDC.mp3', T_sim=T_sim, offset=0.3)
# Ambigutity function plot
fig = plot_fm_ambiguity(fm_sig, fs=fs, frf=100, max_speed=1000, max_range=100, dyn_range=plot_dyn_range)
fig.update_layout(
autosize=False,
width=fig_w,
height=fig_h,
title={'text': 'FM signal Ambiguity Function, ACDC',
'y':0.9,
'x':0.5,
'xanchor': 'center',
'yanchor': 'top'}
)
# Generate FM signal
fm_sig, fs = generate_fm_signal_from_sound('music_samples/speech.mp3', T_sim=T_sim, offset=0.3)
# Ambigutity function plot
fig = plot_fm_ambiguity(fm_sig, fs=fs, frf=100, max_speed=1000, max_range=100, dyn_range=plot_dyn_range)
fig.update_layout(
autosize=False,
width=fig_w,
height=fig_h,
title={'text': 'FM signal Ambiguity Function, Speech',
'y':0.9,
'x':0.5,
'xanchor': 'center',
'yanchor': 'top'}
)
The above prepared functions are available in the pyapril package as well from version 1.7 [2].
from pyapril.sim import FMSimulator
# Generate FM signal
fm_sig, fs = FMSimulator.generate_fm_signal_from_sound('music_samples/speech.mp3', T_sim=T_sim, offset=0.3)
# Ambigutity function plot
fig = plot_fm_ambiguity(fm_sig, fs=fs, frf=100, max_speed=1000, max_range=100, dyn_range=plot_dyn_range)
fig.update_layout(
autosize=False,
width=fig_w,
height=fig_h,
title={'text': 'FM signal Ambiguity Function, Speech',
'y':0.9,
'x':0.5,
'xanchor': 'center',
'yanchor': 'top'}
)